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ABSTRACT 

Numerical investigation of flow and heat transfer in a rectangular duct with a built-in 
circular tube in presence of delta-winglet type vortex generators has been carried out for 
moderate Reynolds numbers and three different angles of attack of the winglets (24°, 30°, 
34°). The investigation is necessitated by the need to enhance heat transfer in fin-tube 
heat exchangers through generation of streamwise longitudinal vortices. The air cooled 
condenser in geothermal power plants also use fin-tube heat exchangers with circular tubes. 
The size of these heat exchangers is huge and often the cost of the condenser is more than 
one-third of the plant cost. Enhancement of heat transfer on the fin surface could reduce 
the size of the condenser. The strategy for heat transfer enhancement involves introduc- 
tion of strong swirling motion in the flow field. The swirling motion can be generated by 
longitudinal vortices. In the present study, the longitudinal vortices have been created by 
delta winglet type vortex generators with common-flow-down configuration. The vortex 
generators are mounted behind the tubes. An element of heat exchanger has been con- 
sidered for detailed study of the flow structure and heat transfer analysis. A significant 
separation delay and removal of zone of poor heat transfer from near-wake of the tubes, 
culminate as a consequence of employing these vortex generators. The heat transfer in the 
near wake region can be enhanced as high as 230%. Results also show a marked increase in 
overall heat transfer in the channel. The enhancement shows a great promise in reducing 
size of the condensers for geothermal power plants. 
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1 INTRODUCTION 


1.1 Problem Description 

Accomplishment of more efficient heat transfer on the air-side of the air-cooled condensers 
used in geothermal power plant applications, improves the overall performance of all these 
plants. The high heat transfer rate must be supported by a modest additional pressure 
drop across the heat exchangers. One of the efficient means of enhancing heat transfer to 
achieve high performance of the geothermal air-cooled condensers has been investigated 
in this work. 

/ 

In gas-liquid crossflow heat exchangers, fin-tubes are commonly used. In such an 
arrangement, the gas generally flows across the tubes and the liquid flows inside the 
tubes. The fins act as extended surfaces providing the bulk of the heat transfer surface 
area on the air side. Even with the extended surfaces, the dominant thermal resistance 
is on the air side, since condensation (phase-change heat transfer) takes place inside the 
tubes. Therefore, in order to achieve significant heat transfer enhancement, strategies 
must be developed, which result in increased heat transfer coefficients on the fin and the 
tube outer surfaces without a large increase in pressure drop. 

Figure 1.1 shows a schematic diagram of the core region of a fin-tube heat exchanger. 
Numerical investigations in this field have been carried out by Biswas et al. (1994), where 
they have tried to achieve the enhancement of heat transfer from the fln surfaces by 
inducing longitudinal streamwise vortices in the flow flled. These vortices were generated 
with the help of delta- winglet type vortex generators on the flat surface (Figure 1.2). 
The longitudinal vortices are developed along the side edge of the delta-winglets due to 
pressure difference between the front surface facing the flow and the back surface. Such 
vortices are also called streamwise vortices because they have axes collinear with the flow 
direction, and once they get developed, they interact with an otherwise a two-dimensional 
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Liquid 



Figure 1.1 Schematic diagram of core region 
of a fin-tube heat exchanger 



Figurre 1.2 Vortex generators on flat fins 




boundary layer and produce a three-dimensional swirling flow that mixes near wall fluid 
with the free stream. This is a type of kinematic mixing and strongly enhances the 
entrainment of fluid from periphery to the core region of the flow field. It results in 
disruption of thermal boundary layer which causes enhancement of heat transfer The 
additional pressure losses happen to be modest because the form drag for such winglet- 
type slender bodies is low. 

In the present work numerical simulations have been carried out for flow and heat 
transfer in fin-tube crossflow heat exchangers with delta winglets mounted behind the tube 
in common- flow- down arrangement as proposed by Pauley et al. (1988). The grid-mesh 
used is of body-fitting type employing finite volume algorithm of Eswaran and Prakash 
(1998). Figure 1.2 presents a sectional view of the suggested arrangement for fin-tube heat 
exchangers with delta winglets placed behind in common- flow- down configuration. 

1.2 The Scope of Present Work 

In order to achieve the objective outlined above, a numerical investigation has been carried 
out. The full Navier-Stokes equations together with the energy equation are solved in a 
rectangular channel with a built-in circular tube having delta winglets placed in common- 
flow-down configuration in the down stream. A detailed analysis of the flow structure 
together with heat transfer characteristics in such a module are studied. 

1.3 Outline of the Thesis 

In Chapter-1 the genesis of the problem has been discussed. Chapter-2 of the thesis pro- 
vides a review of the literature relevant to understanding of the basic mechanism involved 
in augmentation of heat transfer and the algorithms related to solving the Navier-Stokes 
equations. A concise review of modeling aspects of transient flow has also been focussed 
in Chapter-2. The mathematical formulation of the problem for simulation is presented 
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in Chapter-3. In this chapter, the flow geometry, the governing equations, boundary con- 
ditions and grid generation techniques are discussed. Chapter-4 discusses the solution 
algorithm and all the basic description of the finite volume formulation in detail. Chapter 
5 discusses the results for laminar flow with a detailed study of flow physics and associ- 
ated heat transfer. Chapter-6 includes the concluding remarks and the scope for further 
research. 
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2 LITERATURE REVIEW 


2.1 Introduction 

Considerable amount of experimental, analytical and computational research has been 
carried out on enhancement of heat transfer. The enhancement of heat transfer through 
manipulation of surface geometry has been analyzed by many researchers and practition- 
ers and the concerned literature is presented in this chapter. In order to analyze the flow 
structure and heat transfer in such applications, a detailed computational study is needed. 
In such applications, flow regime may possess either laminar or turbulent characteristics. 
Here important numerical techniques developed in the recent past for solving the conserva- 
tion equations for complex flows and heat transfer have been focussed. Numerical solutions 
for laminar flows and heat transfer have been well developed by various researchers but 
the corresponding investigations for turbulent flows involve more rigorous considerations. 
Present literature survey helps in suggesting the work to be carried out to achieve the 
objectives enumerated in Chapter-1. The review of literature follows in two sections, the 
first presenting an overview of work done by various researchers in the area of enhancement 
of heat transfer. The second section presents a review of investigations pertaining to the 
numerical methods with regard to the task of computing flow fields in complex geometries 
thereby it deals with available scheme for solving complete Navier-Stokes equations. 

2.2 Enhancement of Heat Transfer 

The primary interest of present work is to enhance the heat transfer in the gas side of 
fin-tube heat exchangers (with flat fins). With that as aim, we would like to study different 
investigations related to heat transfer enhancement suitable for the applications proposed 
above. 

Investigation of heat transfer enhancement in channel flows where the rate of heat 


5 



transfer between the fluid and channel walls keeps diminishing due to growth of boundary 
layer as the flow tends to get fully developed, is of special interest. Surface protrusions 
on the channel walls destabilize hydrodynamic boundary layer and disrupt the thermal 
boundary layer and thereby enhance the heat transfer between the flowing fluid and chan- 
nel walls. Applications of such flow configurations are the heat transfer between the gas 
and the fin in the case of gas-liquid fin-tube cross flow heat exchangers and the heat trans- 
fer between the flowing fluid and plates in the case of fin-plate heat exchangers. The use of 
a multi-layered surface geometry for plates can disturb evolution towards a fully developed 
flow. A detailed performance data for louvered fin surfaces is provided by Achaichia and 
Cowell (1988). In using louvered fins, enhancement in heat transfer is obtained at the 
price of high-pressure drop and so to take care of that, protrusions in the form of slender 
delta-wings or winglets can be mounted on the fin surfaces (Figure 1.2). As shown, the 
base of the wing remains attached to the fin and the apex faces the incoming stream with 
an angle of attack in this configuration. 

2.2.1 Vortex Generators 

Detailed investigations of the enhancement potential of vortex generators have been carried 
out in the last two decades ( Fiebig et al. (1986), Turk and Junkhan (1986), Eibeck and 
Eaton (1986), Fiebig et al, (1989), Amon (1989), Dong (1989), Zhang (1989), Tigglebeck 
(1992), Yanagihara and Torii (1990), Fiebig et al. (1990), Fiebig et al. (1991), Sanchez 
(1989), Valencia (1992), Biswas et al. (1994) and Valencia (1996)). The idea of using 
vortex generators in compact heat exchangers is being used by the industry. Transverse 
vortex generators, such as ribs and corrugations in channels generate vortices whose axes 
are mainly transverse to the primary flow direction, whereas, longitudinal vortex genera- 
tors , such as delta wings and rectangular wings generate vortices having axes mainly along 
the primary flow direction. The experimental and numerical investigations conducted so 
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far show that longitudinal vortex generators are preferable to transverse vortex generators 
for improvement of heat transfer surfaces in the heat exchangers when both heat transfer 
enhancement and flow loss are taken into account. 

The generation of longitudinal vortices along the side edge of the wing-shaped vortex 
generators appears due to the pressure difference between the front surface facing the flow 
and the back surface. In a channel, these vortices can be made to disturb the growth of 
thermal boundary layer by mixing fluid from the near wall region and the channel core 
region. They can serve to enhance the heat transfer rate while producing less increase in 
pressure loss penalty.The enhancement in heat transfer is caused by the mixing induced 
by these vortices, of the mid-stream and the boundary layer fluids which causes steeper 
temperature gradients to develop near the wall. The application of longitudinal vortices 
for boundary layer control is well known (Pearcey 1961) and such vortex generators find 
prominent use in aircraft industry. 

The study of stream wise longitudinal vortices behind a slender aerodynamic object is a 
research topic of much interest for many years. Hummel and Srinivasan (1967) have made 
important contribution in revealing the complex flow structure behind a delta wing. They 
have conducted experiments and presented pressure distribution and vortex structure of 
the flow around a delta wing of unit aspect ratio with an angle of attack of 20°. Fiebig et al. 
(1986), Fiebig et al. (1991) and Tiggelbeck et al. (1992) have also carried out experimental 
investigations for augmentation of heat transfer by means of longitudinal vortices. Fiebig 
et al. (1986) in their experimental study compared the performance of different kinds 
of vortex generators viz. delta wing, rectangular wing,delta winglet pair and rectangular 
winglet pair in the Reynolds number of 1360 and 2270. Their observation depicts that 
the delta wing is the best vortex generator from the heat transfer point of view. Another 
important feature of their observation is that the heat transfer coefficient increases with 
the increase in angle of attack till the vortex break down takes place. Tigglebeck et al. 
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(1992) used multiple rows of vortex generators in an aligned arrangement within a channel 
and observed their influence on flow structure. They found that the flow structure in the 
wake of the second row has similar qualitative characteristics to that of the first row. They 
also observed that the peak value of the spanwise averaged Nusselt number at the location 
of the wake of second row is strongly dependent on the spacing between the two rows. 
Fiebig et al. (1989) and Biswas and Chattopadhyay (1992) carried out computational 
studies for geometry of delta wing placed inside a channel. They discussed the influence 
of angle of attack and Reynolds number on velocity and temperature fields in the laminar 
range. 

The gas side heat transfer coefficient in a gas liquid fin-tube heat exchanger is small 
as compared to that of liquid side. The mechanism of heat transfer between the gas and 
the solid surface in such cases is to be understood in detail. Consider flow past a tube in 
a channel with neighboring fins. The emerging flow field consists of a horseshoe vortex 
system, a dead water zone at the junction of the tube and the plate and a von-Karman 
vortex street in the middle. To enhance heat transfer in such flow configuration, vortex 
generators can be mounted on the plates, which make the flow field extremely complex. 
Experiments conducted by Dong (1989) and Valenica et al. (1996) show the influence of 
winglet type vortex generators in a channel with a built-in circular tube. Biswaa et al. 
(1994) have observed that in the absence of winglet type vortex generators, relatively little 
heat transfer takes place in the downstream of the circular tube which is the recirculation 
region with low velocity. However, they observed an enhancement of heat transfer as 
high as 240 percent in the wake region behind the cylinder in the presence of winglet 
type longitudinal vortex generators. Fluid flow and heat transfer characteristics through 
a parallel plate channel with transverse ribs, which act as vortex generators, have been 
investigated by Webb and Ramadhyani (1985). Significant heat transfer augmentation is 
obtained in their study. Acharya et al. (1991) explained vortex interaction between two 
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such ribs in a channel for the sub -harmonic excitation of shear layers. Even positioning 
of vortex generators at key locations in the flow has its influence on enhancement of heat 
transfer (Myrum et al. 1992). 

The vortex generators are special surfaces that are used to generate secondary flow or 
vortices by swirling and destabilizing the flow. In addition, these vortices have a tendency 
to last for a long period and to have small cores, allowing an influence on heat transfer 
over long and narrow regions (Fiebig 1998). The vortex generators can generate trans- 
verse or longitudinal vortices depending on the geometry and angle of attack of the vortex 
generators. The axis of transverse vortices (also called spanwise or primary vortices) is 
perpendicular to the main flow direction. These type of vortices are generated from flow 
separation and unstable shear layers, and the flow can remain two-dimensional. The Kar- 
man vortex street in the wake of bluff bodies, e.g., circular and rectangular cylinders, is 
one known example of transverse vortices as studied by Sohankar et al. (1999) and So- 
hankar et al. (2000). Sohankar and Davidson (2001) have numerically investigated the 
effect of inclined vortex generators on heat transfer enhancement in a three-dimensional 
channel. They have confirmed that the downwash of the secondary flow induced by vortex 
generators increases mixing of cold and hot fluid and enhances heat transfer. They also 
showed that when thicker vortex generators are used stronger and bigger streamwise vor- 
tices are formed downstream of the vortex generators which give rise to a higher Nusselt 
number. Fiebig (1998) concluded that in general, longitudinal vortices are more efficient 
than transverse vortices for heat transfer enhancement. In addition, longitudinal vortices 
enhance locally and globally in steady flow, while transverse vortices generate a negligible 
global heat transfer enhancement in steady flow. 

Numerical investigation of convective heat transfer and pressure drop in wavy ducts is 
carried out by Blomerius and Mitra (2000) in the laminar to transitional range of Reynolds 
number (600-2000). They have determined the geometrical parameters for the best ratio 


9 



of heat transfer to pressure drop for three-dimensional channels consisting of corrugated 
plates. Both 2-D and 3-D geometries are considered by them where two-dimensional com- 
putations show the onset of self-sustained flow oscillation above a critical Reynolds number 
and an associated enhancement of heat transfer.Two-dimensional parametric studies have 
been used to identify potentially interesting three-dimensional geometries. 

2.3 Numerical Methods for Solving Navier-Stokes Equation 

For last three decades various solution techniques have been proposed for full Navier- 
Stokes equations. The solution techniques for incompressible flows do not have any ex- 
plicit equation for pressure and so due to spatial coupling of pressure and velocity it 
becomes quite difficult to obtain direct numerical solutions. Pressure does not have the 
usual thermodynamical meaning in incompressible flow problems and becomes a relative 
variable, which adjusts itself instantaneously for the condition of zero mass divergence 
to be satisfied at all computational cells. The primitive variable formulation of the in- 
compressible Navier-Stokes equations has following most distinctive feature. Because in 
incompressible Navier-Stokes equations speed of sound is treated to be infinite, the pres- 
sure field cannot be calculated by explicit time advancement procedure and so it require at 
least a partially implicit determination which takes into account the coupling between the 
pressure and velocity fields as well as the effects of the velocity boundary conditions. This 
pressure field can be conveniently computed using stream function-vorticity approach for 
two-dimensional flows but it again creates problem for computation of three-dimensional 
flows due to absence of single scalar stream function. 

Primitive variable approach can be employed for three-dimensional problems which 
can be followed without encountering non-physical wiggles in pressure distribution. Har- 
low and Welch (1965) have used a staggered grid for the dependent variables in their well 
known MAC (Marker and Cell) method which is one of the earliest and widely used explicit 
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methods for solving the full Navier-Stokes equations. In this method, solution of velocities 
are obtained using a two step procedure. In the first step called the predictor step, the 
provisional values of velocity components are computed explicitly using advection, diffu- 
sion and pressure gradient of the earlier time steps. This explicitly advanced provisional 
velocity field may not necessarily ensure a divergence free velocity field. Hence, in the sec- 
ond step called the corrector step, the pressure and velocity components are corrected so 
that the velocity field satisfies the continuity equation. This is done through the solution 
of a Poisson equation for pressure. A relative technique, known as pseudo-compressibility 
method, developed by Chorin (1967) involves a simultaneous iteration of the pressure and 
velocity components. Vicelli (1971) has shown that the method due to Chorin (1967) and 
the MAC method are equivalent. 

Since implicit methods have no such restrictions, they are more attractive. Patankar 
and Spalding (1972) have introduced an efficient method known as SIMPLE (Semi-Implicit 
method for Pressure Linked Equations). This method is based on a finite- volume dis- 
cretization of the governing equations on a staggered grid. In order to improve the con- 
vergence involved in the pressure- velocity coupling, several variants of SIMPLE algorithm 
have been developed. The SIMPLER algorithm of Patankar (1981) and the SIMPLEC 
algorithm of Van Doormaal and Raithby (1984) are improvements on SIMPLE. Although 
the changes to incorporate SIMPLEC into SIMPLE algorithm are minor, the consequences 
can be as great as it eliminates the approximations made in SIMPLE while deriving the 
pressure-velocity corrections. 

Application of finite volume methods using non-orthogonal coordinates and collocated 
grids is reported by Rhie and Chow (1983) and Peric (1985). Ferziger and Peric (1999) have 
shown that collocated arrangement converges faster than the staggered variable arrange- 
ment and has advantages when extensions such as multigrid techniques and non-orthogonal 
grids are considered. Mukhopadyay et al. (1993) has developed a numerical method for 
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incompressible geometries. Integral mass and momentum conservation equations are dis- 
cretized into algebraic form through numerical quadrature. The physical domain is divided 
into a number of non-orthogonal control volumes which are iso-parametrically mapped on 
to standard rectangular cells. Numerical integration for unsteady momentum equation is 
performed over such non-orthogonal cells. The results exhibit good accuracy and justify 
the applicability of the algorithm. Kobayashi and Pereira (1991) have modified the momen- 
tum interpolation method suggested by Peric (1985) and named it as Pressure- Weighted 
Interpolation Method-corrected (PWIMC). In this method, non-orthogonal terms in the 
momentum equations were solved explicitly, whereas in the pressure-corrections they were 
dropped. The finite volume algorithm developed by Eswaran and Prakash (1998) for 
solving time dependent three-dimensional full Navier-Stokes equation has the following 
advantages: 

• It uses collocated grid arrangement, hence all the advantages of collocated grid 
arrangement are gainfully utilized. 

• The essence of the algorithm is SMAC so it is suitable for modeling unsteady flows 
(see Kim and Benson, 1992). 

• The governing equations are discretized in the physical plane itself without coordi- 
nate transformation. 

This algorithm has been used to solve the flow problem and extended to solve the energy 
equation for the present work. Peric et al. (1988) carried out a comparative study of finite 
volume numerical methods with staggered and collocated grids. They concluded that the 
collocated method possesses no disadvantage relative to the staggered grid version. The 
treatment of boundary conditions and the implementation of higher order differencing 
schemes is simpler. In some cases, the collocated version provides faster convergence. It 
behaves similar to the staggered version under the variation of grid and under-relaxation 
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parameters, and provides solutions of equal accuracy at the same or lower computing cost. 
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3 FORMULATION OF THE PROBLEM 


3.1 Introduction 

The need to augment the heat transfer in fin-tube heat exchangers so as to reduce their 
size and thereby cost, has led us to compute the flow and heat transfer characteristics in 
a horizontal channel with a built-in circular tube and a pair of winglets. This numeri- 
cal study can provide information about the heat transfer enhancement in fin-tube heat 
exchangers and effect of various parameters on the transport enhancement. 

3.2 Statement of the Problem 

The computational domain is shown in Figure 3.1. Two neighboring fins form a channel 
of height H, width B = 11.25H and length L = 25H. The circular tube of diameter D = 
2.8H is located at a distance Li = 6.25D from the inlet. The blockage ratio, D/B has been 
kept 0.25. The position of the winglet is = 7.5H and Xb = 10.5H, Ya = 2.8H and Yb 
= 1.125H, h = 0.67H and angle of attack has been varied. Air has been used as working 
fluid, hence the Prandtl number of the study is 0.7. 

The flow field in the channel with built-in circular tube and a winglet pair is charac- 
terized by the following parameters: 

• Reynolds number 

• Blockage ratio (D/B) 

• Position of the cylinder in the channel (Li/L) 

• Velocity profile at the channel inlet 

• Channel height (H/B) 

• Angle of attack of the winglets (P) 
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Figure 31 Computational Domain 



• Height of the winglets (h/H) 


• Position of the winglets (X^/L, Y^/L) 

3.3 Governing Equations 

For laminar flow of an arbitrary spatial control volume V bound by a closed surface S 
the three-dimensional Navier-Stokes equation can be expressed in the following general 
convection-diffusion-source integral form: 


dt 



pu 


■d't 


(3.1) 


dt 



(plt<p 


r^V4>) -dt = f s^dv 

Jv 


(3.2) 


Where p represents the fluid density, u is the fluid velocity, 4> stands for any vector 
component or scalar quantity, S^is the volumetric source term. For incompressible flow of 
Newtonian fluid, the equation becomes 


dt 


j) it ■d't = Q 
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S^dV 
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(3.4) 


and the source term for momentum equation becomes f pl.d~^ where I is the unit 
tensor. 

Here we work with Cartesian components of velocity. So <f) can be the three-Cartesian 
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components of velocity u, v, and w as well as any scalar e.g., temperature which need to 
be determined. It is to be noted that Equation (3.2) is a perfectly general equation from 
which the three fundamental equations of flow and heat transfer (i.e., continuity, momen- 
tum and energy) can be derived. The Table 3.1 displays substitutions for (j) in general 
transport equation to obtain the conservation equations for mass, momentum and energy. 
Basu et al. (2001) have computed flow and heat transfer in a cross flow past a circular 
tube placed in a channel. 


Table 3.1: Values of the variables in the general transport Equation (3.2). 


Equation 

<f> 



Continuity 

1 

0 

0 

Momentum 

Uj 


dp 

dxj 

Energy 

T 

k 

Cp 

0 


3.4 Boundary Conditions 

The governing differential equations are elliptic in space and parabolic in time. The 
boundary conditions for all confining surfaces are required which are outlined below. 

• Top and bottom plates 

u = V = w — 0 (No slip boundary condition) and = 0 
T = Ty, (Tm represents wall temperature) 


• Side Walls 

|^ = ^ = 0,u = 0 (Free slip boundary condition) and || = 0 
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• Channel inlet 


u = Uoo, V = w = 0 and 1^ = 0 


T = Too 


• Channel exit 

^ + Uav^ = 0 (Orlanski boundary condition) 

(where (j) represents either of u, v, w or T) 

P = Poo 

• Obstacles (Cylinder and Winglet) 

u = V = w = 0 and ^ = 0 (where n signifies the normal direction) 
T = T^ 
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4 GRID GENERATION AND METHOD OF SOLU- 


TION 

4.1 Introduction 

The numerical solution of the governing equations in the prescribed computational do- 
main with given initial and boundary conditions involves two major discretization steps. 
First being the discretization of the domain with the help of a numerical grid used for 
computation. The governing equations are discretized on the basis of this numerical grid. 
Various ways of discretization can be applied to the above set of partial differential equa- 
tions. Among these the most common methods are those using finite differences, finite 
elements and finite volumes. In the present work a finite volume scheme is chosen for spa- 
tial discretization of the governing equations. As a result, a conservative formulation of 
the discretization is possible, assuring automatically the satisfaction of the global balance 
of the conserved quantities, independent of the coarseness of the numerical grid. 

4.2 The Grid 

Figure 4.1 shows the schematic of the two dimensional grid used in the present work. 
The grid generation technique uses block partitioning method to generate the grid. The 
whole domain is divided into two blocks. In first block where winglets have been placed, 

cartesian grid is generated. In the second block where a built-in circular tube is placed, 

/ 

a differential equation technique is used for generation of the body-fitting grid. The two 
dimesional grid-mesh is stacked in the third direction to obtain the three dimensional 
grid-mesh. The grid generation technique adopted for body-fitting type is described in 
brief in the following sections. Three categories of physical domains may exist viz. (1) a 
simply-connected domain, (2) a doubly-connected domain and (3) a multiply-connected 
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Figure 4.1 The Grid System 





domain. Heie simply-connected domain is used in the body- fitting region. 

4.2.1 Grid Generation Technique 

The two-dimensional grid has been generated using tranfinite interpolation and differ- 
ential equation method. The method of tranfinite interpolation essentially uses linear 
interpolation scheme to compute the interior points by using the coordinate values from 
the boundaries. The resulting algebraic grid is not smooth in general because the slope 
discontinuities at the boundary propagate in the interior of the domain. That is why the 
mesh obtaiiUKl by algebraic mapping is further improved by the use of Partial Differential 
Equations (PDE) technique in which a system of PDE is solved for obtaining the location 
of the grid points in the physical space. But in the computational space the transformed 
grid is of rectangular shape with uniform spacing. We can use any of the three types of 
partial differential equations, namely elliptic, parabolic and hyperbolic where elliptic equa- 
tions are preferred for closed geometries, hyperbolic equations are used for domains with 
outer boundary not prescribed and similarly parabolic equations are used where boundary 
is closed on one side (prescribed) and open on the other. The present grid is generated 
using of elliptic partial differential equations. In any grid generation technique geomet- 
ric information is fetched from the boundaries. Hence the steady state boundary valued 
nature of elliptic equations makes them most favorite and so they are widely used grid 
generation. A Laplace equation or a Poission equation with a Dirchlet boundary condition 
can be used for this purpose. Standard Poission equation is of the form 


v2e = p 


V^r] = Q 


(4.1) 
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Where P and Q are known as control functions. 

If P = Q = 0 then, the equations get modified to Laplace form as 

V2^ = 0 


V^77 = 0 (4.2) 

The above equations can be solved by finite difference technique to get the location of 
the interior points from the boundary. The initial guess needed to solve the Laplace 
equation can be generated using algebraic mapping. This Laplacian operator can provide 
quite smooth grids. The grid lines remain equally spaced if boundary curvature is absent 
but near curved boundaries they show a tendency to concentrate. This drawback can be 
controlled with appropriate selection of control functions P and Q. Also control functions 
can be objectively employed for enforcing orthogonality over the selected edges. That 
may simplify implementation of boundary conditions, but use of control functions remains 
associated with additional complexity in the transformed domain. 

4.2.2 Computational Procedure for Grid Generation 

The algorithm used for the generation of the present grid can be outlined as follows: 

1. Input is fed in the form of geometrical description. 

2. The computational grid is defined on the curvilinear ^ - r] plane based on the number 
of grid points. 

3. Boundary points along the edges of the physical domain are defined. 

4. Now an initial guess is generated for the Laplacian by algebraic grid generation 
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method. This step uses just a linear interpolation scheme to compute the interior 
points. 

5. The guess is improved a little by adopting simple four-point formula. 

6. The actual Laplacian is in the form of 

- 2bx^r, + cxnn = 0 


"h cyjff) — 0 


(4.3) 


where 



h = x^Xj^ + y^yrj 


is solved till convergence is obtained. Here subscripts indicate partial derivatives of 
Cartesian variables with respect to curvilinear variables. 

The truncation error in numerical solutions is very much dependent on the nature of grid 
because it is here that the descretized approximation of differential equations is solved. 
\Nq target for such a grid which gives minimum truncation error, and that requires testing 
of certain properties of the grid which have strong effect on the error. The parameters 
that affect these properties can be used to express the quality of the grid. These are 

• Transformation Jacobian 

• Skewness 
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• Aspect ratio 


• Adjacent cell ratio 

A desiiable giid should have nonzero tranformation Jacobian almost at each location of 
the giid cells (they should be almost of similar order of magnitude with no quite small 
values appealing anywhere), least possible skewness, aspect ratio close to 6, and adjacent 
cell ratio close to one. 

4.3 Finite Volume Method 

The governing conservation equations are discretized using the finite volume approach 
of Eswaran and Prakash (1998). Here the solution domain is divided into a number of 
contiguous (finite) control volumes (CV) which are defined by the coordinates of their 
vertices assumed to be connected by straight lines. The coordinates of the control volume 
vertices are calculated by grid generation procedure. A collocated grid arrangement is 
employed and all the dependent variables(u, v, w, p, T) are defined at the same location 
- the centroid of the control volume (Figure 4.2). E, W, N, S, T and B indicate the 
six neighboring control volume centers for the east, west, north, south, top and bottom 
neighbors. The face center points e, w, n, s t and b are located at the intersection of the 
lines joining the mid points of opposite edges. For example, te, be, ne. se are the midpoints 
of the edges that form the east face with e as the center of the cell face. The advantages 
of collocated grids over staggered grids have already been pointed out in Chapter-2 as 
explained by Peric et al. (1988). Before the time stepping can be done it is necessary to 
calculate geometrical parameters which include the surface vectors for the six faces e, w, 
n, s, t and b for each finite volume and its volume. 
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4.4 Surfaces and Volumes 


The finite volume vertices are numbered 1 to 8 in the manner shown in Figure 4.2. The 
outward surface normals and volume can be found in the following manner as suggested 
by Kordulla and Vinokur (1983) and Eswaran et. al.(1995). 

Defining = ry - Vj where n and Vj are the position vectors of points i and j respec- 
tively, we have expressions for surface area vectors of various faces given as follows 

• East face 


= ^(■^74 X T^ss) 

• West face 

= ^(T^i6 X T^52) 

• North face 

= ^(■^27 X T^63) 


• South face 



(4.4) 


(4.5) 


(4.6) 


(4.7) 


• Top face 



(4.8) 
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• Bottom face 


^6 = ^(T^13XT>24) (4.9) 

The volume of any cell is calculated using the cell coordinates, with the assumption 
that the linear segments to form the six cell faces join the cell corners (Eswaran and 
Prakash, 1998).The expression for volume of any cell becomes 

= + (4.10) 

4.5 Discretization Procedure 

The main steps of the discretization procedure to calculate the convection and the diffusion 
fluxes and source terms are outlined below. The rates of change of integral fluxes and the 
source terms are integrated over the cell volume, whereas the convection and diffusion 
terms form the sum of fluxes through the faces of the control volume. 

4.5.1 Discretization of Continuity Equation 
Equation (3.1) is discretized in the following way: 

j pit -dt ^ p{^ =Y P'^j ■ 

^ j=e,w,n,s,t,b ^ 

where S/is the surface vector representing the area of the j th cell face and is the velocity 
defined at the face center j. 

The discretized form of the continuity equation becomes 

YFj = Fe + F^ + Fn + Fs + Ft + Fi, = 0 (4.12) 

j 
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where Fjis the outward mass flux through face j, defined by 



4.5.2 Discretization of General Equation 
(a) Rate of Change 


The value of the dependent variable (j) at the centroid of the control volume (the geo- 
metric center) represents an average over the CV as a whole. Thus 


dt 



{p(i>V)T - ip4>v)l 

At 


pV 


At 


( 4 . 13 ) 


where V is the volume of the cell. 


(b) Convection Fluxes 

The surface integral over convection flux of variable (j) can be approximated in the following 
form 




where </>jis the value of 4> at the center of face j. Thus 


( 4 . 14 ) 


/ P^lt • - 1 - FyjCj)^ + Fn4)n + Fg^g + Ft4>t + Fb4>b ( 4 - 15 ) 

Js 

where 4>e is the interpolated value of the variable cj) at the east face center etc. This can be 
evaluated by using a central difference linear interpolation between the neighboring nodal 
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values (j}p and (pp. At east face the value of <f>e, is given by 


(j^e — 




(j)p + 


Vp 

Vp + Vp 


4>b 


(4.16) 


where V'b and Vp are volumes of the cells around the points E and P respectively and 
and (j>p are the values of the dependent variables at these points. In a collocated 
grid system, all dependent variables are defined at the same location hence exactly the 
same interpolation scheme is used to express all of them at the interfaces. The central 
difference approximation to compute the convection flux may lead to numerical stability 
problems. Therefore the convection flux is blended with a first order upwind differencing 
scheme (UDS), and the difference between the central difference scheme (CDS) and UDS 
approximation as 


F,4 = (F.W''“ + 7|(-F'.'^.) 


\CDS 


{Fe4>e) 


UDS 


(4.17) 


The upwind differencing scheme is based on assumption that the convected cell face is 
equal to that at the upstream cell along the same coordinate direction. Thus the value 4>e 
at the east face is assigned the value (j)p ifue > 0, i.e., the flux is negative. This can be 
conveniently summarized as 

= [[P’..0]]-fe[[-F.,0]l 

The symbol [[ , ]] signifies greater of the two quantities enclosed inside the brackets. 
Similar expressions can be written for the other cell faces (Eswaran et al. 1995). 
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Fw4^w — 4^p 0 ]] 4^e 0 ]] 


+T'{^” + -'^f[[^"».oj) + fe((-f^^ ( 4 . 19 ) 

F„<i>n = i'p ((^;, 0]] - <t,K ([-F„ 0 )] 

FA. = 't>p [[i=’..0]l-?is[|-J’„0)] 

(vs + Vp'^^ Vs + Vp^^} ” [[“-^e'0]]| (4-21) 

i='«4i = 0p[[Ft,Oll-.^p[[-.F„Oll 

.F10» = (6pl[n.O]]-fe[l-.Fj,Ol] 

In a fully implicit method the upwind parts are implicit and they are incorporated in 
the coefficients of the unknown velocity during the pressure-velocity iteration. The CDS 
terms on the other hand are evaluated using the previous iteration values and used as a 
source term on the right side of the same equation. This is the so called deferred correction 
approach of Khosla and Rubin (1974). Multiplication of the explicit part of 7 (0 < 7 < 1 ) 
allows the introduction of numerical diffusion. For the first order upwind difference scheme 
7 = 0 and for the second order central difference 7 = 1 . The deferred correction approach 
enhances the diagonal dominance of the coefficient matrix. However, the present solution 
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scheme being explicit, the accuracy depends on the value of 7 . 


(c) Diffusion Fluxes 

Peric (1985) developed the expression for the diffusion fluxes in the flnite volume approach 
for first time. However, the exact expression in three dimensions are not explicitly available 
in open literature. Eswaran and Prakash (1998) evaluate the diffusion flux of a variable <j) 
through the cell faces in the following manner 





For any face we can write 


(4.24) 


j — aiin} H- a2n^ + (4.25) 

.<<*s 

where and are any three linearly independent (not necessarily orthogonal) unit 

vectors. Therefore 


V</> 


y,-=v^-( 


axn} + 0:2 


(4.26) 


= aiV<j> ■ + az^(j> ■■n^ 

If A0® are the differences in ^ between the two ends of the line segments AX\ 

AX^, AX3 then 

A<f>^ = V(f> ■ A^, A<^2 ^ • A^ (4.27) 

If A^, A^, A^ are in the directions tP, and respectively then it follows from 
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Equation (4.27) that 


__ 


/\({> 


A(^3 




(4.28) 


where AA\ AX^&nd AX^ are the magnitudes of AA^, A^and A^.Combining Equa- 
tions (4.26) and (4.28) we have 




■■ ai 




+ 0i2- 


Axi 


• + 0:3 


AX2 


(4.29) 


To obtain ai, ors and 0:3 we express 


v} = (nn,ni2,ni3) 


n2 = (n2i,n22,n23) (4.30) 

✓N. 

V? = (n3i,n32,n33) 

where nn, ni 2 , riiz are the Cartesian components of and which can be easily determined 
where AXl^ AX], AX] are the three components of vector A^. The 
other values n2\ ... n^z can be similarly obtained. Hence equation (4.22) can be written as 
a matrix equation involving (i, j = 1,3) forming a known coefficient matrix of unknown 
column vector of 0 : 1 , o; 2 and 0:3 related to Cartesian components of the surface vector j 
as S\j , , cLTidi Szj . 

Using Cramers rule in such a matrix equation we can solve for ai, a^and az as follows 
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(4.31) 


Di 

D’ 




where D is the determinant of the coefficient matrix. Di is obtained by replacing the first 
column of D by the column with the elements 5y, S^j and Szj. 

The diffusion flux is made of two distinct parts - normal derivative diffusion flux and 
cross derivative diffusion flux. The second part arises from the non-orthogonality of the 
grid. The normal derivative diffusion flux is treated implicitly and is coupled with the 
implicit part of the convective flux to calculate the main coefficients of the discretized 
equations while the cross derivative diffusion flux is treated explicitly to avoid the possi- 
bility of producing negative coefficients in an implicit treatment. This term together with 
explicit part of convective flux is added to the source term. The example of east face is 
taken to illustrate the diffusion model (Figure 4.3). Given the edge center values (pte, (t>be, 
<l>sei <Pne WO can get the normal diffusion term and the cross diffusion term 

and Finally the diffusion flux is computed by 



ar 


I ^i)e . 0ne 

AJ¥> * AX2 


(4.32) 


In order to find the edge center values for evaluation of cross derivative diffusion flux, 
the following interpolation scheme is proposed. 


4>te = 


Vtb 1 , ^p a _l 

YJ—VP + TT-(PTE + 
Vtot Vtot 


Vr 

w 


(l>E + 


tot 


Ve 

Vtot 


4>t 


(4.33) 


Vtot = VtE + Vp -I- Vx + V£; 
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Figure 4.3 Face representation to illustrate 
diffusion model 



(4.34) 


, Vbe , , Vp , , , I j 

4>be — Tr~9P + JJ—9BE + T7~'9b + JT-CPB 
not not not not 


Vtot — ^BE + ^P + + Ve 


ytot ytot Vtot ytot 


(4.35) 


Vtot = Vne + Vp + Vn + Ve 


j, _ ^SE , . ^P 1 I A \ A 

9se = Tr~9P + T7—9SE + T7~9e + Tr~9S 
Vtot vtot Vtot Vtot 


(4.36) 


Vtot = VsE -\-Vp -\-Vs + Ve 


where is the volume of the cell that is located at the top east of the cell P. The other 
edge center values of the dependent variables can be interpolated in a similar way. 

(d) Sources 

The source term is to be integrated over the cell volume. The contribution of a source 
term in the discretized equation may be thought of 
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(4.37) 


S^dV ^ {S^)p V 

Apart from the real source Stp, explicitly treated parts of the convection and diffusion fluxes 
maybe added to S^. The momentum equations contain pressure terms. These terms are 
also treated explicitly. Its discretization is analogous to that of the ordinary diffusion flux, 
i.e., for the i th momentum equation the pressure term is 



— / pfii • d~^ fti ^ PjSij (4.38) 

j 

where pj is the pressure at the j th face center and Sijis the i th component of the source 
vector for face j. 

4.6 Pressure and Velocity Coupling 

The velocity and pressure field satisfying the mass and momentum conservation laws can 
be found using a procedure akin to Simplified Marker and Cell (SMAC) method. This 
method offers an efficient and easy way of pressure velocity coupling. It is basically a 
semi explicit method. The momentum equations are discretized in an explicit manner 
with the exception of the pressure gradient terms that are treated implicitly, and the 
continuity equations are also enforced implicitly. This can be expressed mathematically 
by the following two discretized equations of momentum and continuity 

+ + = (4-39) 

3 j 
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and 


= 0 (4.40) 

j 

The momentum equations are solved using the guess values of velocity and pressure 
field. The provisional velocity components are calculated from the following equation 

where u°p is the value of velocity at earlier iteration and Su is the pressure term. This 
provisional velocity in general will not satisfy the continuity equation. The continuity 
Equation (4.40) in another form reads as follows: 


E “ 

3 3 

/where Fj is the uncorrected mass flux obtained from the provisional velocities and Fj is 
the mass flux correction. To evaluate the terms on the right side of the Equation (4.42) it is 
necessary that the variables (velocities and pressure) are to be known at the cell faces. Due 
to the non-staggered arrangement if the variables at the cell faces are evaluated by linear 
interpolation between the adjacent cell center quantities then the pressure iteration does 
not converge. The solution leads to a Checker Board pressure field (i.e., spurious oscillation 
of pressure may occur) as shown by Rhie (1981). This problem can be circumvented by 
using the concept of Momentum Interpolation (Majumdar, 1988). The essence of the 
concept is that the velocities at the cell faces are computed by linear interpolation of 
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the convective and diffusive terms but not the pressure term. Thus the method with the 
collocated variable arrangement relies indirectly on the staggering idea. Following this 
idea the interpolated velocity at the east face of the control volume is obtained in the 
following manner 


Ue = (up , ve) - 



(4.43) 


where 


Vp = U°p 


At 

Wp 


(n + Ft) 


(4.44) 


= + Fi) (4.45) 

and the over bar indicates a linear interpolation using equation (4.16). So the uncorrected 
mass flux for the east face using Equation (4.43) becomes 

F: = pit, ■ te = p{ltp , its) • te - AtVp • te - (4.46) 

Equation (4.46) in its generalized form for any cell face is in the for 


F* = pltj ■ 'fj - AtVpj ■ (4.47) 

To enforce mass conservation, the velocity and pressure corrections are introduced, linked 
by 



(4.48) 
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This step simulates the effect of a staggered grid (face-center velocity and cell-center 
pressure) and reaps the benefits of faster convergence. Prom Equation (4.48) we get the 
mass flux correction at the face j 

Ej = pu'j • Sj = -AtWj ■ S'j (4.49) 

Since Equation (4.49) is exactly of the same form as equation for diffusion flux Equation 
(4.24), with only the variable interchanged, therefore the same method can be used for 
computing Fy Combining Equations (4.42) (4.47), and (4.49) yields 

—AtVp'j ■ j = pV^j ■ j — AtVpj ■ 

or 

AtVpJ • = 0 (4.50) 

(5 „ = ptj ■ y, - Atvp, . y,) 

which gives the pressure correction. After the solution of the pressure correction equation, 
the nodal velocity, mass fluxes and pressure are updated. The corrected pressure is 
obtained by 


pu+l ^pn+p' 


(4.51) 


Instead of solving Equation (4.39) again with the corrected pressure field, velocity calcu- 
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lations can be corrected in the fore 


(4-52) 

The corrected velocity field is obtained by 

(4.53) 

Using this new velocity distribution, the energy equation is discretized and solved. On 
integrating the conservative form of the steady state energy equation over a control volume 
(neglecting the dissipative terms as it is significant for high speed flows), we arrive at the 
following temperature equation 

a (VT . y ) = (pV, . y,) T, (4.54) 

i ^ j 

which is solved iteratively to get the temperature distribution of the current time step. 

The first terms of Equation (4.50) and Equation (4.54) are expanded following the 
philosophy that was followed during the discretization of Equation (4.24). Subsequent 
regrouping of the normal derivative terms yields the standard Poisson Equation. The 
Poisson equation can be solved by a multitude of methods. We have used Gauss Seidel 
iteration method. This method is quite reliable and robust. The flow chart with the main 
steps of the algorithm is given in Figure 4.4. 

4.7 Solution Algorithm 

The velocity, pressure and temperature are calculated bj' the following procedure of 
Eswaran and Prakash (1998). The procedure can be summarized in the following way: 

1. Make initial guess of the pressure, velocity and temperature fields. Use Equations 
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(4.41) and (4.44) to calculate cell center velocities up and vp. use linear interpolation 
{Equation(4.16)} to obtain the face center quantities Vj. 

2. Compute the mass flux through cell face j using equation (4.47) i.e., 

f ; = pitj ■ t, - Atvp, ■ 

3. Use equation (4.48) to compute the flux correction at the cell face j using Equation 
(4.49) i.e., 


= -AtVp; • 


This is computed using the formulation for diffusion fluxes [Equations (4.24) - (4.36)]. 
For calculating Fj , the variable (p is to be replaced by p'. 

4. Compute the residue, 3? for each cell, 

5. Calculate the cell-center pressure correction (using Gauss-Seidel method) from the 
relation 


Pp=Pp-\-UJ 


ttpAt 


where to is the relaxation factor and Cp stands for the diagonal coefficients of Equation 
(4.50) , which is calculated from 


ai 

Oil 

“2 

Q!2 

0=3 

as n 

LAXi 

VJ AX^ 

e AX2 

n AX2 

5 AX3 

CO 

<1 

1 
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where ai, etc. are the same as that are used in Equations (4.24 - 4.36). 

6. If > e go to step 3. 

7. Store the updated mass flux through the cell faces using Equation (4.47) 

F* = ■ 'tj - AtVpj ■ 'fj 

8. Store the updated pressure at the center of the cells 

Pp = Pp+ p'p 

9. Store the cell center corrected velocity 

Up = Up-\- Up 

10. Calculate the residual for each cell 

3 3 

11. Calculate the cell center temperature using Gauss-Seidel method, to obtain 

Tp = Tp + fir 

bpa 

where Q is the relaxation factor and bp stands for the diagonal coefficients of Equa- 
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tion (4.54) i.e., calculated from the relation 


bp = —ap 


12. If 3firms > e go to 10. 

13. Go to step 1 and repeat the process until the steady state is reached. It can be 
shown that solving Equation (4.42) is same as solving 



At 


and finding the corrections 





such that 


= u: + 


Thus the procedure above can be called solution of Pressure Linked Equations. 
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Figure 4.4 Iterative Solution Scheme 










4.8 Numerical Stability Considerations 

The Semi-Explicit method used here relies on explicit differences and suffers from time 
step restrictions. For a given mesh the choice of the time step is determined through 
stability analysis which has to take care of two conditions. First, fluid should not be 
allowed to cross more than one cell in any one time step. This restriction is derived from 
the Courant-Friedrich-Lewy (CFL) criterion of stability given by 


5ti < min 


f5x 

Sy 

5z 1 

1 w ’ 

|z;i’ 

w j 


(4.55) 


where the minimum is with respect to every cell in the domain. Typically 5t is chosen 
equal to one-third to two-third of the minimum cell transient time. 

Secondly, when a nonzero value of kinematic viscosity is used, momentum must not 
diffuse more than one cell in one time step. A linear stability analysis shows that the 
restrictions on grid Fourier number in this case will yield 


1 {6xf {5yf {5z) 

2 {5xf + {5yf + (Szf 


(4.56) 


Finally the minimum of the above two time increments is chosen for computation. 
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5 RESULTS AND DISCUSSION 


5.1 Introduction 

A three-dimensional numerical solution of flow and heat transfer for a plane channel with 
built-in circular tube and delta winglets in common flow down arrangement is presented 
in this section. A 51 x 41 x 17 grid-mesh is used and the computations have been carried 
out for Reynolds number of 500, 1000 and 1410 and at different angles of attack. The 
divergence free criterion is satisfied using an upper bound of 10"'^. The blockage ratio 
(D/B) is fixed at 0.25. Air has been taken as the working fluid, hence the Prandtl number 
of the present study is 0.7. The heat transfer characteristics have been compared on the 
basis of spanwise averaged Nusselt number ( Nus ) based on bulk mean temperature(Tf,). 
The entire analysis in this section is based on the time-averaged flow and temperature 
fields. 

5.2 Flow and Heat Transfer Characteristics in a Rectangular Chan- 
nel with built-in Circular Tube 

In this section detailed analysis of flow and heat transfer characteristics for rectangular 
channel with a built-in circular tube is presented. The results presented here are for a 
Reynolds number of 1000. 

5.2.1 Flow Characteristics 

Figure 5.1 shows the streamline plot on the horizontal mid plane of the channel i.e. at 
z = H/2. The figure clearly shows the two symmetrical standing vortices at location A. 
These are the von Karman vortices. The flow separation in the wake region of the tube is 
clearly established. 

The limiting streamline plot in the region close to the plate is shown in Figure 5.2. A 
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Figure 5. 1 Streamline plot on horizontal midplane 



Figure 5.2: Limiting Streamlines on a plane close to the bottom plate 



saddle point of separation and a horseshoe vortex system is observed from the figure. The 
incoming flow does not separate in the traditional sense but reaches a stagnation or saddle 
point of separation (marked A) and goes around the body. The nodal point of attachment 
(marked C) and the separation lines which form circular arcs across the tube are seen in 
the figure. The flow above the lower wall hits the front of the tube. A significant part of 
it moves downward and creates a region of reversed flow in front of the stagnation line. 
On each side of the tube, one finds a region of converging streamlines (marked G). These 
are the traces of horse shoe vortices. Behind the body one finds two areas of swirling flow 
(marked E) which are the footprints of an arch vortex. There is a wake stagnation point 
(marked F) further downstream of the body. 

Figure 5.3 shows streamline plot on the vertical midplane. The wake of the cylinder 
exhibits a strong three-dimensional behavior. The strong normal velocity component, w is 
caused by pressure gradient in the vertical direction in the wake region. The figure clearly 
shows the footprints of the horseshoe vortices (marked as Si and S 2 ) and the reattachment 
line or the wake stagnation line (marked S 3 ). 

The horseshoe vortices (Si and S 2 ) are initiated near the junction of cylinder and 
the channel walls. The formation of these vortices is displayed in Figure 5.4. As the 
fluid approaches the stagnation line of the cylinder, its velocity head contributes to total 
pressure. The velocity in the boundary layer on the channel walls increases in the vertical 
direction away from walls. This results in stagnation pressure gradient. This causes flow 
towards the walls which interacts with the main stream. The fluid rolls up forming vortices 
which are swept around the cylinder base and are carried down stream. This results in 
formation of horseshoe vortices which form a pair of counter rotating longitudinal vortices. 


49 




Cylinder 


Plate 



Karman vortex 
street 


Horseshoe 
vortex system 


Dead water zone 


Figure 5.4 Flow around a cylinder on a flat plate 
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5.2.2 Heat Transfer Characteristics 

Figure 5.5 shows the distribution of spanwise average Nusselt number based on bulk mean 
temperature (i^) for rectangular channel and rectangular channel with buiH-in circular 
tube. The Nug at the leading edge, for both cases, is very high due to largo temperature 
difference between the incoming colder fluid and the hot channel walls. rectangular 

channel Nug decreases continuously in downstream direction. The velocity u da y 

1 • , . , , ,.-11 Mpnce, same fluid 

layer is predominantly two-dimensional, and the cross flow is negligible, 

particles are in contact with the channel walls. As a result, the temperature 

, ^ flow direction, 

between the fluid and channel walls decreases and hence Nug decreases m 

_. r icine channel with 

Figure 5.6 shows iso-Nusselt number plot on the bottom plate tor piaii 

1 - 1 . . . 1 ..f Nusselt number 

built-in circular tube. Figures 5.5 and 5.6 show an abrupt increase or i'' 

, , n vortex system 

at the location of circular tube. This is due to formation of horsesnoe 

, . ^ in Figure 5.4. 

that consists of two counter-rotating longitudinal vortices as already snowi 

_ . , i„vpr and produce 

These vortices interact with predominantly two-dimensional boundary 

, , 1 iUz, frpe stream. This 

a three-dimensional swirling flow that mixes the near wall fluid with tne 

f Vipat transfer and 

leads to disruption of thermal boundary layer causing enhancement oi ne 

+11 he is a separated 

increase in span averaged Nusselt number. The wake zone behind the tuu 

, f transfer in this 

dead water zone. It has fluid recirculating at very low velocity and so neai/ 
zone is quite poor. 

5.3 Flow Field in Plane Channel with a Built-in Circnlnr Tube 
and Winglet Pair 

1 a built-in 

This section presents a detailed analysis of flow in a rectangular channe 

1 j number of 1000. 

circular tube and winglet pair. The results correspond to a Reynoia& 

• +q 1 midpl^^® the 

Figure 5.7 shows the streamline plot of the flow field on the horizontal 
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Re = 1000, Pr = 0.7 
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Figure 5.5: Nug for Plane Channel and Channel with Circular Tube 




Figure 5.6: Iso-Nusselt number plot on the bottom plate 





Figure 5.7: Streamlines on the horizontal mid-plane of the channel 


channel at z = Hf2. The flow separation and the wake formation are clearly visible. The 
wake is symmetric at this stage. Point A shows the saddle point of separation and Wi,W 2 
show the location of winglets. Figure 5.8 shows the limiting streamlines for the velocity 
fleld on a horizontal plane close to the bottom wall for the Reynolds number of 1000. The 
winglet pair brings about a swirling motion, which has a dominant component of transverse 
motion. The transverse momentum transfer to otherwise separated boundary layer delays 
the separation. A saddle point of separation and a horseshoe vortex system are observed 
in the figure. The incoming flow reaches a stagnation point or a saddle point of separation 
(marked A) and goes around the body. The nodal point of attachment (marked C) and 
separation lines are also visible. On each side of the tube, one finds a region of converging 
streamlines (marked G). These are the traces of the horseshoe vortices. Behind the tube 
one finds two zones of swirling flow (marked E) which are foot prints of an arch vortex. 
The wake stagnation point F is located at a distance slightly more than one diameter 
in the downstream. The winglets are located at Wi and W 2 . The twisted streamlines, 
indicating the swirling motion are visible around the winglets. 

5.3.1 Vortex Formation by Delta Winglets 

Figure 5.9 shows the vortices generated by the winglet. It generates two types of vortices, 
namely, main vortices and horseshoe vortices. Horseshoe vortices are induced due to 
velocity variation in boundary layer over the base plate. The mechanism is similar to 
the formation of the horseshoe vortex system generated in flow past a circular tube as 
explained in section 5.2.1. Main vortices are created due to pressure difference between 
the pressure side and the suction side of the winglets. 

Figure 5.10 (a) , (b), (c) show cross-stream velocity vectors at various axial locations in 
the channel with the circular tube and the winglet pair as the built-in obstacles. Schematic 
location of the winglet pair is shown in Figure 3.1. The view-window of this figure is near 


56 




close to the bottom plate 





Figure 5.10 Streamline plot in cross plane at various X locations 

(a)X=ll, (b)X=16, (c)X = 21 







the exit plane. The axial location of Figure 5.10(a) is at the immediate downstream of the 
trailing edge of the winglet pair. The velocity vectors depict the secondary flow pattern 
on various cross-stream planes in the channel. The results correspond to a Reynolds 
number of 1000 and /? = 30°. The pressure difference between pressure surface and suction 
surface of the winglet induces a complex vortex system. The vortex system generally 
behaves like a free vortex motion (apparent from the velocity vectors) , with the exception 
of the core region where a rotational (forced) vortex motion persists. The vortices are 
also longitudinal vortices, since the axes of the vortices are parallel to the main flow 
direction. The longitudinal vortices are centered at Mi and M 2 . The swirl, imposed by 
the longitudinal vortices, serves in creation of downwash on the bottom plate. The figure 
clearly reveals that the secondary flow on the cross stream planes describes a common-flow- 
down situation. Mendez et al. (1998) have explained the kinematics of flow and mechanism 
of transport enhancement due to the longitudinal vortices. The strength of the vortex 
motion decreases in the downstream due to viscous resistance initiated at the confining 
walls. Vortex interaction with common flow between the vortices directed towards the wall 
has been extensively studied by Pauley and Eaton(1998). The downwash causes strongest 
distortion of the boundary layer over the greatest streamwise distance. Figure 5.11 displays 
the basic flow structure for the swirling nature of the flow in both common-flow-up and 
common- flow- down configurations. The resulting downwash fountaion in common-flow- 
down configuration may be observed in the figure for our case of counter-rotating vortex 
pair. 

Figures 5.12 (a), (b) and (c) represent the streamlines on the cross stream plane at 
different axial locations, namely at X=ll, X=16 and X=21 in the channel. Here Mi and 
M 2 are the main vortices generated by the winglets. Also are shown the horseshoe vortices 
Hi and H 2 generated by winglets. These vortices are also called corner vortices (Biswas 
et al. (1996)). The streamlines confirm the existence of a strong vortical motion. The 
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Downwash 



Figure 5.11 Flow structure due to delta winglets 



Figure 5.12 Streamline plot in cross plane at various X locations 

(a)X = ll, (b)X=16, (c)X = 21 


centers of main vortices are located at a distance of ±0.35 from the vertical mid plane. 
Induced vortices become prominent as the pair of longitudinal vortices travel in the down- 
stream 

5,4 Heat Transfer Characteristics in Plane Channel with Built-in 
Circular Tube and Winglet Pair 

Figure 5.13 compares span-averaged Nusselt number distribution on the bottom wall in a 
channel for the case of built-in circular tube, with and without winglet pair. The Nusselt 
number distribution starts with a high value and decreases in the flow direction due to the 
growth of boundary layer on the bottom plate. The Nusselt number assumes a high value 
near the axial location of the leading edge of the circular tube. As the fluid approaches the 
stagnation line of the circular tube, it slows down and its pressure increases. The smaller 
velocity within the boundary layer in the vicinity of the bottom plate, which supports the 
circular tube, leads to a smaller pressure increase. Thus the induced pressure gradient 
causes a downward flow towards the bottom wall. Having impinged on the plate, the 
fluid rolls up forming vortices which finally wrap around the front half of the tube and 
extends to the rear of the tube (Goldstein and Karni(1984)). The resulting motion can be 
best described as a horseshoe vortex system. The horseshoe vortices are also longitudinal 
vortices and they promote disruption of thermal boundary layer. The Nusselt number is 
enhanced due to the horseshoe vortices. 

An abrupt increase of Nusselt number is observed near the axial location of trailing 
edges of the winglet pair. Two counter-rotating longitudinal vortices culminate during 
creation of a complex vortex system due to the winglet pair. The formation of longitudinal 
vortices brings about better mixing and the Nusselt number becomes high near the local^ion 
of the winglets. The span-averaged Nusselt number reaches a high value (9.31) near the 
axial location of the leading edge of the circular tube. This is primarily governed by the 
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Common-flow-down configuration 
Re = 1000, Pr = 0.7, P = 30° 
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Figure 5.13 NUg for plane channel, channel with circular tube 
and channel with circular tube and winglet pair 



formation of the horseshoe vortices at the tube and bottom wall junction as explained 
earlier. Up to the location of the leading edge of the winglet pair, the distribution of span- 
averaged Nusselt number is same for the cases with and without winglet pair. However, 
beyond that location span-averaged Nusselt number increases for the case of channel with 
the tube and winglet pair. The peak value of span averaged Nusselt number takes place in 
the neighborhood of the trailing edge (x=10.5) of the winglet pair and its value is 10.74. 
The maximum enhancement in heat transfer (230%) is also observed at the same location. 
At the exit of the channel, the span-averaged Nusselt number is 6.03. The span-averaged 
Nusselt number at the exit of the plane channel is 3.26. The average Nusselt number value 
of the entire bottom plate for the case of a channel with built-in circular tube and the 
winglet pair is 40 percent more than that for the case of a channel with a built-in circular 
tube. 

Iso-Nusselt number plot on bottom plate is shown in Figure 5.14. From entry till the 
end of the circular tube, the variation of Nu is similar to that in Figure 5.6. In Figure 5.14 
winglets are placed at locations Wi and W 2 . Figure shows increase in Nusselt number 
just after the location of winglet pair. 

Figure 5.15 shows the distribution of local Nusselt number along the width of the 
channel at an axial location of x=ll. This location is in the immediate downstream of 
the trailing edge of the winglet pair. The flow field and the temperature field used for 
calculation of local Nusselt number are at any arbitrary time instant. The peaks ( Nu 
= 19.2) reconcile the effect of counter rotating longitudinal vortex pair. Maxima of the 
local Nusselt number occur at the spanwise locations where the downwash due to vortical 
motions impinges on the bottom plate. 


65 




o 


Figure 5.14 Iso-Nusselt number plot on the bottom plate 
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5.5 Effect of Reynolds number and Angle of Attack on Heat 
Transfer 

Figure 5.16 compares the variation of span averaged Nusselt number for a Reynolds number 
1000 and at different angles of attack. The span-averaged Nusselt number ceases to be a 
strong function of the angle of attack of the winglet pair beyond /9 = 24°. 

Variation of Nu^ for various Reynolds numbers at same angle of attack is shown in 
figure 5.17. It shows that Nus increases with Reynolds number. For plane channel, it is 
known that Nug increases with increase in Re. So even for the case of circular tube with 
winglet pair the correlations may still show increase of Nug with Re (as it is observed 
here) but their exact nature needs further investigation. 

5.6 Model Validation 

The model validation was performed through comparison with the experimental results 
of Eckert and Soehngen (1952) and numerical results of Chun and Boehm (1989) for a 
Reynolds number of 200. Figure 5.18 shows the comparison of local Nusselt numbers 
along the circumference of the tube. Present computation compares favorably with the 
experiments of Eckert and Soehngen (1952). 
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Common-flow-down configuration 
Re = 1000, Pr = 0.7 



Figure 5.16 Variation of NUs with the angle of attack ( (3 ) 






Re = 200 ; Pr= 0.7 
Eckert and Soehngen (1 952) 
Chun and Boehm (1 989) 
Present Study (at the midplane) 



Figure 5. 18 Comparision of local Nusselt number distribution on the Tube Surface 



6 CONCLUSION AND SCOPE FOR FUTURE WORK 


6.1 Conclusion 

A three dimensional numerical study on the flow and heat transfer characteristics in plane 
channel with built-in circular tube and winglet pair in cross flow of the tube has been con- 
ducted. The delta winglet pair has been placed here in common flow down configuration. 
The duct was designed to simulate a passage formed by two neighboring fins in a fin-tube 
heat exchanger. The governing conservation equations of fluid flow are derived in integral 
form. An explicit control volume formulation devised by Eswaran and Prakash (1998) is 
used to solve the equations. A body fitting grid has been used for simulation of flow and 
heat transfer and investigation has been carried out for full domain. In computations of 
Biswas et al. (1994) the consideration of half domain for computation imposes a symmetry 
about vertical mid plane. This leads to suppression of possible vortex shedding. In present 
case of full domain computation, any chance of vortex shedding fully persists. 

The streamline plot near the bottom wall clearly indicates the three-dimensionality 
of flow field, where all nodal points of attachment and saddle points of separation can 
be seen. The strength of longitudinal vortices beyond delta winglet pair is observed in 
cross plane and is seen to decrease in downstream direction. The consequence of this 
variation of vortical strength is reflected in decrease of span-averaged Nusselt number. 
The distribution of Nu and iso-Nusselt number shows the high heat transfer zone behind 
the tube and beyond the trailing edge of winglet pair. This is the poor heat transfer zone 
due to presence of dead water wake zone behind the tube (Basu et al. 2000) in absence of 
winglet pair. 

Delta winglet type vortex generators, placed in the aft of the circular tube, cause 
separation delay. This brings about narrowing of the wake and induces a strong swirling 
motion. The characteristic feature of the swirling motion is appearance of counter-rotating 
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vortical structures on the cross-stream planes behind the circular tube. The counter- 
rotating motion takes the cooler fluid from the underside of the core of the vortices, swirls 
around to upper side, entraining warmer fluid from the outside into the core of the vortices. 
As a result, the core region gets heated up and at the same time the heat transfer from the 
wall to the fluid is enhanced due to disruption of thermal boundary layer. Thus the zone 
of poor heat transfer from the near wake of the circular tube is removed. The creation of 
streamwise longitudinal vortices entail 230 percent enhancement of span-averaged Nusselt 
number in the near wake. The augmentation of average Nusselt number on the bottom 
plate is about 40 percent. 

6.2 Scope for Future Work 

The results of this work validate the unsteady flow modes as a useful technique in improv- 
ing heat transfer efficiency for air-cooled heat exchangers. Similar to the common flow 
down conflguration analyzed here the case of common flow up configuration could also be 
taken up and a comparison made with present results. The latter case will have leading 
edges of winglets more separated than trailing edges. One thing is obviously expected that 
latter case will reduce the wake zone significantly. But the behavior of Nu variation etc. 
is still a question. Also the circular tube could be replaced by oval tube and delta winglet 
pair to be mounted in the aft of the tube and results compared with present case. Even 
more than a pair of delta winglets could be mounted. In all these case delta winglet pair 
could well be replaced by rectangular winglet pair and the results compared. 

Also the variation of vortical strength with downstream distance beyond winglet pair 
could be investigated and its correlation with Nu could be confirmed. Establishing opti- 
mum angle of attack for these cases is still a potential problem to be tackled. 
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